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ABSTRACT 


Multi-dimensional  upwind-differencing  schemes  for  the  Euler  equations  are  reviewed.  On 
the  basis  of  the  first-order  upwind  scheme  for  a  one-dimensional  convection  equation  the  two 
approaches  to  upwind  differencing  are  discussed:  the  fluctuation  approach  and  the  finite- 
volume  approach.  The  usual  extension  of  the  finite-volume  method  to  the  multi-dimensional 
Euler  equations  is  not  entirely  satisfactory,  because  the  direction  of  wave  propagation  is 
always  assumed  to  be  normal  to  the  cell  faces.  This  leads  to  smearing  of  shock  and  shear 
waves  when  these  are  not  grid-aligned.  Multi-directional  methods,  in  which  upwind-biased 
fluxes  are  computed  in  a  frame  aligned  with  a  dominant  wave,  overcome  this  problem,  but 
at  the  expense  of  robustness.  The  same  is  true  for  the  schemes  incorporating  a  multi¬ 
dimensional  wave  model  not  based  on  multi-dimensional  data  but  on  an  ’’educated  guess” 
of  what  they  could  be. 

The  fluctuation  approach  offers  the  best  possibilities  for  the  development  of  genuinely 
multi-dimensional  upwind  schemes.  Three  building  blocks  are  needed  for  such  schemes: 
a  wave  model,  a  way  to  achieve  conservation,  and  a  compact  convection  scheme.  Recent 
advances  in  each  of  these  components  are  discussed;  putting  them  all  together  is  the  present 
focus  of  a  worldwide  research  effort.  Some  numerical  results  are  presented,  illustrating  the 
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1  Introduction 


CFD  algorithms  for  the  coming  generation  of  massively  parallel  computers  will  have 
to  be  extremely  robust.  They  will  most  likely  be  implemented  on  adaptive  unstruc¬ 
tured  grids,  and  will  be  used  for  ambitious  simulations  of  steady  and  unsteady  three- 
dimensional  flows.  In  such  a  complex  environment  there  is  little  place  left  for  hand¬ 
tuning  parameters  that  regulate  accuracy,  stability  and  convergence  of  the  computa¬ 
tions.  A  typical  algorithm  will  make  very  intensive  use  of  local  data,  with  a  minimum 
of  message  passing. 

Algorithms  of  this  nature  exist  already  in  CFD:  they  are  the  upwind-differencing 
schemes,  computationally  intensive  but  unsurpassed  in  their  combination  of  accuracy 
and  robustness.  While  these  favorable  properties  are  explainable  for  one-dimensional 
methods,  it  is  a  stroke  of  luck  that  upwind  schemes  work  as  well  as  they  do  for 
two-  and  three-dimensional  flow.  Their  design  is  commonly  based  on  one-dimensional 
])hysics,  namely,  the  solution  of  the  one-dimensional  Riemann  problem  that  describes 
the  interaction  of  two  fluid  cells  by  finite-amplitude  waves  moving  normal  to  their 
interface.  The  inadequacy  of  this  technique  clearly  shows  up  when  the  numerical 
.solution  contains  shock  or  shear  waves  not  aligned  with  the  grid,  for  instance,  by  a 
loss  of  resolution. 

The  need  to  incorporate  genuinely  multi-dimensional  physics  in  upwind  algorithms 
was  recognized  as  early  as  198.'5  by  Phil  Roe  [1].  A  study  of  discrete  multi-dimensional 
wave  models  by  Roe  followed  in  1985  (ICASE  Report  85-18,  also  [2]).  but  it  took  until 
1991  [3]  before  any  algorithms  based  on  such  wave  models  became  truly  successful. 
Important  contributions  to  this  development  were  made  by  Herman  Deconinck  and 
collaiforators  [3,  4]  at  the  Von  Karman  Institute  in  Brussels.  The  new  upwind  schemes 
are  formulated  on  unstructured  grids  with  data  in  the  vertices  of  triangular  or  tetra¬ 
hedral  cells. 

While  genuinely  imdti-dimensional  methods  were  slowly  developing,  partial  suc- 
ces.ses  were  booked  by  putting  some  multi-dimensional  information  into  the  Riemann 
solvers  used  in  conventional  upwind  schemes.  In  particular,  it  became  the  fashion 
to  obtain  a  plausible  wave-propagation  angle  from  the  data,  rather  than  accepting 
the  angle  dictated  by  the  grid  geometry.  The  earliest  work  of  this  kind  is  due  to 
.Steve  Davis  [5];  it  recently  was  picked  up  by  a  number  of  authors:  Levy,  Pow’ell  and 
Van  Leer  [6],  [7],  Dadone  and  CIrossman  [8,  9],  Obayashi  and  Coorjian  [10],  Tamura 
and  F'lijii  [11].  Roughly  speaking,  they  apply  Riemann  solvers  in  several,  physically 
appealing,  directions;  I  shall  refer  to  their  work  as  the  multi- directional  approach. 

Related,  but  closer  to  the  genuinely  multi-dimensional  approach  is  the  work  of 
Rumsey,  Van  Leer  and  Roe  [12.  13,  14.  15]  and  Parpiaand  Michalek  [16,  17].  These  au¬ 
thors  indepemlently  developed  almost  identical  multi-<limensional  wave  models  based 
on  minimizing  wave  strengths.  These  wave  models  requires  only  two  input  states,  just 
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Figure  1;  Two  views  on  scalar  upwiiul  ditfereTiciiig:  (a)  nodal-point  interpretation; 
(1))  finite- volume  interpretation. 

as  a  regular  Riemann  solver. 

In  support  of  these  quasi-multi-dimensional  approaches,  aimed  at  putting  better 
physics  into  interface  fluxes,  some  authors  have  dedicated  efforts  to  improving  the 
interpolation  or  reconstruction  step  that  precedes  the  flux  calculation.  On  a  struc¬ 
tured  grid  the  reconstruction  of  a  non-oscillatory  distribution  of  flow  variables  from 
their  cell-averages  usually  is  done  dimension  by  dimension:  a  fully  multi-dimensional 
reconst rucion  is  indispensable  in  achieving  higher  accuracy.  Barth  and  Frederickson 
'IS  indicated  how  to  reconstruct  a  smooth  function  up  to  arbitrarily  high  order  on 
an  unstructured  triangulation;  .Abgrall  [19]  showed  how  to  implement  truly  multi¬ 
dimensional  limiting  of  higher  derivatives. 

In  this  lecture  I  shall  review  a  decade  of  efforts  toward  multi-dimensional  upwind- 
differencing.  with  the  accent  on  the  very  latest  developments.  The  discussion  is  limited 
to  the  multi-dimensional  physics  that  goes  into  these  methods:  multi-dimensional 
reconstruction  will  not  be  further  mentioned.  For  a  somewhat  different  emphasis 
or  point  of  view  the  reader  is  referred  to  three  excellent  other  reviews  of  multi¬ 
dimensional  methods  [’JO.  21.  4]  that  have  been  presented  in  the  past  year. 


2  Two  views  of  one-dimensional  upwinding 

In  order  to  appreciate  the  problems  surrounding  multi-dimensional  unwinding  it  is 
useful  to  consider  the  principles  of  one-dimensional  upwinding.  The  reader  is  assumed 
to  be  familiar  with  the  theory  of  conservative  upwind  schemes:  as  a  tutorial  Roe's  [22] 
review  article  is  recommended. 

r|)wind  differencing  is  a  way  of  differencing  convection  terms.  For  the  scalar 
coiuection  equation 

Ut  +  cui  =  0.  ( 1 ) 

the  simplest  upwiiul-difference  scheme,  of  first-ordei  accuracy,  reads 
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Figure  2:  Two  approaches  to  upwind  differencing  for  the  Euler  equations:  (a)  fluctu¬ 
ation  approach:  (b)  finite-volume  approach. 
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■Scheme  (2,. 3)  can  be  regarded  as  a  formula  for  updating,  from  to  either  the 
nodal-point  value  of  u  in  x,,  or  the  cell  average  of  u  in  cell  i.  These  two  view-points 
are  illustrated  in  Figures  la  and  lb.  The  distinction  is  significant,  because  it  leads 
to  distinct  methods  for  more  complex  equations.  In  the  development  of  schemes  for 
the  one-dimensional  Euler  equations,  the  first  view-point  has  led  to  the  concept  of 
fluctuation  splitting,  due  to  Roe  [23.  22];  the  second  view-point  is  that  of  Godunov 
[24]  and  has  led  to  the  projection/evolution  or  reconstruction/ evolution  concept  of 
finite- volume  schemes,  due  to  Van  Leer  [25.  26.  27].  Below  1  shall  review  the  formulas 
pertinent  to  each  approach. 


2.1  Fluctuation  splitting 

.Assume  the  system 

-P  F((')x  =  0  (4) 

represents  the  Euler  equations  in  conservation  form.  i.e..  U  =  [p.  pu,  pE)^  is  the 
vector  of  conserved  state  quantities  and  F[U)  =  [pu.pu^  4-  p,puH)^  is  the  vector  of 
their  fluxes.  The  equation  shows  that  any  local  imbalance  of  the  fluxes  causes  the 
local  solution  to  change  in  time.  .Such  a  local  imbalance  is  called  a  fluctuation  by  Roe 
[28.  Ij.  If  source  terms  are  present,  their  value  must  be  included  in  the  fluctuation 

[32]. 

Define  the  matrix  c\((/)  as  the  derivative  of  F{U)  with  respect  to  (/,  so  that 


dF{U)  =  A{U)dU.  (5) 

It  is  essential  for  the  technique  of  fluctuation  splitting  that  this  differential  relation 
be  replaced  by  an  exact  finite-difference  analogue,  namely, 

AF  =  AAU,  (6) 

where  A  indicates  a  difference  between  neighboring  nodal  points.  Roe  [29]  has  in¬ 
dicated  how  to  construct  a  mean  value  /I  of  A  such  that  Eq.  (6)  holds  exactly  for 

arl)itrary  pairs  of  state  vectors.  For  a  calorically  perfect  gas  a  suitable  mean  value 
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T 

can  ea.sily  l)e  olitained  l)v  introducing  the  the  parameter  vector  w  =  ^  a.  H)  . 
Since  butli  I  (tc)  and  F{w)  are  quadratic  in  the  components  of  w.  it  follows  that  Eq. 
(6)  is  satisfied  by  A  =  ,4  (f/(?e)),  where  w  is  the  algebraic  average  of  w. 

Fluctuation  splitting  requires  that  the  matrix  .4  be  split  into  its  positive  and 
negative  parts,  i.e.. 

A  =  A-^  +  a-.  (7) 

so  that 

AF  =  .4+Af/  +  A- AC.  (8) 

A  popular  name  for  this  procedure  is  ”flux-difFerence  splitting”;  the  term  ‘‘fluctuation 
splitting"  is  preferable  because  it  includes  source-term  splitting.  The  first  term  on  the 
right-hand  side  combines  disturbances  that  propagate  forward:  in  consequence,  this 
term  is  used  to  update  the  right  nodal  point.  The  second  term  combines  backward- 
moving  disturbances  and  is  used  to  update  the  left  nodal  point.  This  concept  is 
illustrated  in  Figure  2a.  Conservation  is  ensured  because  the  two  terms  add  up  to  a 
perfect  flux  difference.  The  first-order  update  formula  becomes 

- 1?  ■  (9) 

In  practice  it  often  pays  to  abandon  the  matrix  notation  and  expand  AF  and  AF 
in  terms  of  the  individual  disturbances.  This  yields 
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AF  =  '£akRk. 

k=l 

3 

(10) 

AF  =  ^  X^akRk, 

(11) 

where  A;,  is  an  eigenvalue  of  A,  /?*.  is  the  corresponding  eigenvector,  and  is  the 
wave  strength:  note  that  Eqs.  (fi)  and  (10)  imply  Eq.  (11).  By  considering  that  each 
fluctuation  may  move  forward  or  backward  through  the  grid,  we  recover  the  splitting 
formula  (8): 


AF 


^  XhOkflk  +  XkCtkRk 

.\fc<0  .\fc>0 


k=l  k=\ 

A^au  -b  A-au. 


(12) 


2.2  Finite-volume  approach 

In  the  finit<-- volume  approach  the  focus  is  on  the  numerical  flux  function  F{Ui,Ur), 
a  recipe  for  computing  the  interface  fluxes  from  the  states  Fi,  and  Fr'  on  the  left 
and  right  sides  of  the  interface.  The  generic  formula  for  updating  cell  averages  of  the 
conserved  quantities  is 
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(13) 
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In  (lodunov's  first-order  scheme  the  interface  Hux  is  taken  from  the  solution  at 
<  >  0  of  Riemaun’s  initial-value  problem  with  input  data 

U{x,0)  =  ill,  X  >  0, 

O'ix,  0)  =  Uft,  X  <  0; 

this  is  illustrated  in  Figure  2b. 

For  many  applications  it  is  not  necessary  to  use  the  exact  solution  to  this  problem, 
hence  the  activity  in  the  research  area  of  “approximate  Riemaim  solvers”  [23,  30,  31]. 
Adopting  Roe’s  [23]  approximate  solution,  which  is  the  exact  solution  of  the  locally 
linearized  equation 

Ut  +  Au,==0,  (16) 

we  find  three  equivalent  formulas  for  the  interface  flux: 

F{UL,Un)  =  Fl  +  A-AU,  (17) 

F{UL,Un)  =  Fn-A^AU,  (18) 

F{Ul,Ur)  =  ^-{Fl  +  Fr)-\A\AU,  (19) 

where 

=  (20) 

In  practice  the  formula  ( 19)  is  preferred  because  of  its  symmetry;  the  expanded  form 
is  ^ 

F{U^,Ur)  =  Ufl  +  Ffi)  -  i  f:  ]A,]a,F,.  (21) 

^  A;=l 

Inserting  the  flux  (19)  into  the  finite-volume  scheme  (13)  yields  an  scheme  that, 
with  the  help  of  the  identity  (6).  reduces  precisely  to  the  fluctuation  scheme  (9).  Yet. 
there  exists  an  important  difference  between  Eqs.  (19,13)  and  Eq.  (9):  in  the  latter 
the  matrix  A  must  satisfy  the  identity  (6)  in  order  to  maintain  conservation,  while 
in  Eq.  (19)  the  matrix  ]/A]  may  be  derived  from  ariy  average  A  without  endangering 
conservation.  The  flux  formula  (19),  due  to  Van  Leer  [32,  33],  preceded  the  fluctuation 
approach  of  Roe  [23],  based  on  (6),  by  a  decade. 


(14) 

(16) 


3  Intermezzo:  how  good  is  one-dimensional  up- 
winding? 

To  appreciate  the  superior  accuracy  and  robustness  of  upwind  differencing  in  one 
dimension,  consider  the  numerical  results  shown  in  Figure  3  and  4,  taken  from  [34] 
and  [35],  respectively.  In  Figure  3a  the  exact  and  discrete  Mach-number  distributions 
for  choked  flow  through  a  converging-diverging  channel  are  superimposed.  First- 
order  fluctuation  splitting  was  used,  including  source-term  splitting  [22,  36]  and  a 
special  splitting  near  the  sonic  point  [34].  .Although  the  update  formula  is  only  first- 
order  accurate,  it  can  be  shown  that  the  scheme  yields  second-order  accurate  steady 
solutions.  In  fact,  in  the  steady  state  the  scheme  reduces  to  the  two-point  box  scheme 
on  all  meshes  except  near  a  sonic  point  and  inside  a  shock  structure,  where  it  becomes 
a  three-point  scheme.  This  yields  the  smooth  transition  through  the  sonic  point  and 


rlie  crisp  shock  fraiisitiou  in  the  displayetl  results.  I'ie;nre  ;}1)  shows  the  resirlua.1- 
coiiveroeiice  histories  for  tliree  increasin,s,Iy  powertull  marching  techniques;  global 
time-stepping,  local  time-stepping  and  characteristic  time-stepping  [dd];  these  look 
mie\entful.  In  Figure  4a  a  shockless  trausonii'  soliititjn  is  reached  from  initial  \alues 
containing  7  shocks  and  8  sonic  points:  again,  the  residual-convergence  history  in 
Figure  4b  for  local  time-stepping  shows  nothing  unusual. 

It  is  this  type  of  performance  we  wish  to  preserve  when  e.xtending  upwind  differ¬ 
encing  to  higher  dimensions. 


Euler  Equations  for  Channel  Flow  Euler  Equations  for  Channel  Flow 


Figure  .'5;  (i'hoked  How  through  a  converging-diverging  channel.  com])Uted  with  a 
Huctuation  scheme,  (a)  Initial  and  final  .\Iach-number  distributions:  (b)  residual- 
convergence  histories  for  global,  local  and  characteristic  time-stepping. 


Eul«r  Equatjont  for  Channel  Flow  Convergence  History 

Mach  Number  in  Channel 


Figure  4:  Transonic  flow  in  a  converging-diverging  channel,  computed  with  a  fluctua- 
iion  scheme,  (a)  Initial  and  final  .Vlach-number  distrilmtions;  (h)  residual-convergence 
history  for  local  time-stepping. 
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Figure  5:  Extending  the  finite-volume  method  to  two  dimensions  by  solving  one¬ 
dimensional  Riemann  problems  at  all  cell  faces.  The  arrows  symbolize  the  exchange 
of  information  between  cells  in  the  direction  normal  to  their  interface. 

4  Multi-dimensional  extension  of  the  finite- volume 
method 

The  standard  way  to  extend  upwind  differencing  to  the  multi-dimensional  Euler  equa¬ 
tions  is  still  the  same  as  indicated  by  Godunov  et  al.  [37]  in  1961.  For  first-order 
accuracy,  initial  values  are  assumed  to  be  constant  in  each  cell,  just  as  in  one  di¬ 
mension:  fluxes  at  cell  interfaces  again  follow  from  solving  one-dimensional  Riemann 
problems  of  the  type  (14.15).  with  x  now  measuring  distance  along  the  normal  to  the 
interface.  This  is  illustrated  by  Figure  5. 

It  is  the  projection  of  the  true  initial  values  onto  cellwise  constant  distributions 
(or  linear  [25,  26]  or  quadratic  [25,  38,  39]  or  even  higher-order  distributions  [40])  that 
creates  discontinuities  at  the  interfaces.  This  leads  us  to  introducing  plane  wave  fronts 
parallel  to  the  interface,  and  selecting,  out  of  all  possible  directions,  the  interface 
normal  as  the  direction  for  wave  propagation.  If  the  solution  contains  only  shock 
and/or  shear  waves  aligned  or  nearly  aligned  with  the  grid,  this  choice  happens  to  be 
the  correct  one,  and  high  resolution  of  such  waves  can  be  achieved  in  the  steady  state, 
just  as  in  one  dimension.  If,  however,  such  waves  are  far  from  aligned  with  the  grid, 
they  get  misrepresented  by  the  upwind  scheme  as  pairs  of  grid-aligned  waves,  as  shown 
in  Figure  6  for  a  shear  wave.  Thus,  a  grid-oblique  stationary  wave  may  be  represented 
by  several  grid-aligned  running  waves,  leading  to  higher  numerical  dissipation  and  a 
considerable  loss  of  resolution. 

Another  purely  numerical  artifact  caused  by  grid-aligned  upwinding  is  the  pres¬ 
ence  of  pressure  disturbances  across  a  grid-oblique  shear  layer.  First  observed  by 
Venkatakrishnan  [41],  the  explanation  wa.s  provided  by  Rumsey  et  al.  [12];  this  phe¬ 
nomenon  is  further  discussed  in  Section  4.2. 

From  the  above  critique  one  should  not  conclude  that  in  higher  dimensions  the 
standard  upwind  methods  are  inferior  to  other  methods:  the  loss  of  accuracy  just  is 
much  more  obvious  for  upwind  methods. 


R 
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Figure  6:  Misinterpretation  ot  a  grid-oblique  shear  wave  by  grid-aligned  upwinding. 


Figiue  7;  Fluxes  in  a  frame  aligned  with  a  wave  front  oblique  to  the  grid  lines. 

4.1  Multi-directional  methods 

The  smearing  of  oblique  shock  waves  in  numerical  solutions  has  received  considerable 
attention,  and  a  proportionally  large  research  effort  has  been  spent  in  mending  this 
weakness.  The  prevailing  idea  is  to  solve  the  Riemann  problem  in  a  direction  more 
appropriate  than  the  grid  direction.  One  immediate  consequence  of  leaving  the  grid- 
aligned  frame  is  that  solving  one  Riemann  problem  no  longer  suffices.  Figure  7  shows 
that,  in  two  dimensions.  l)oth  flux  vectors  in  the  rotated  frame  are  needed  for  the 
construction  of  the  thixes  normal  tv)  tlu'  inleifacv^. 

Consider,  for  example.  Figure  8,  showing  a  rotated  coordinate  system  aligned  with 
level  lines  representing  a  shock  front  in  a  discrete  solution.  It  makes  sense  to  solve 
a  one-dimensional  Riemann  problem  in  the  direction  normal  to  the  front,  i.e.  using 
the  flow-velocity  components  in  that  direction;  this  yields  the  flux  in  the  normal 
direction.  The  input  states  for  the  Riemann  solver  are  Uix  =  Ui  and  (/rx  =  F'r-  The 
flux  tangential  to  the  shock  should  be  obtained  from  state  values  located  at  L||  and 
R\\\  using  ('i  and  f’/?  once  more  would  completely  destroy  the  effect  of  the  rotation 
[7.  14]).  These  values  could  be  approximated  by 

^  f  "r||  =  -  {l‘L  +  f  r)  :  (22) 

this,  howeviu'.  implies  cf'ntral  dilferencing  along  the  shock  and  leads  to  odd-ev<>n 
decoupling  in  that  direction  [6.  7.  42]. 


S 


Figure  S;  A  simple  multi-directional  flux  formula. 


Figure  9:  Input  states  for  the  Riemann  problems  in  tlie  flux  computation  according 
to  Levy  et  al. 

In  the  work  of  Davis  [5.  4)5] .  dating  back  as  far  as  1983,  the  computation  of  the 
tangent  flux  actually  is  more  complicated  than  that  of  the  normal  flux.  The  more 
recent  work  of  Levy  et  al.  [6,  7.  42]  and  Dadone  and  CIrossman  [8,  9]  is  more 
mature  in  that  the  fluxes  are  treated  without  distinction.  Figure  9  shows  how  pairs 
of  input  states  to  the  two  Riemann  problems,  L’flx )  and  {1'l\\A'r\\)i  are  selected 

according  to  Levy  et  al.  In  their  first-order  method  the  input  states  in  the  rotated 
frame  are  obtained  by  linear  interpolation  between  neighboring  states  in  a  ring  of  cells 
surrounding  the  interface;  Dadone  and  Grossman  simply  take  the  value  in  the  nearest 
cell,  which  apparently  adds  to  the  robustness  of  the  method.  Another,  wider  ring  of 
cells  is  needed  for  achieving  second-order  accuracy. 

Various  choices  can  be  made  for  the  rotation  angle  of  the  frame  in  which  the 
Riemann  problems  are  solved.  A  sensitive  quantity  is  the  direction  of  the  velocity- 
difference  vector.  Vfi  —  Vi,  which  was  adopted  by  Davis  and  also  is  crucial  to  the 
approach  of  Rumsey  and  Parpia  (see  Section  4.2).  Levy  et  al.  use  the  direction  of 
the  velocity-magnitude  gradient  VjV  |.  which  can  detect  both  shock  and  shear  waves, 
while  Dadone  and  Grossman  use  the  pressure  gradient  Vp.  which  only  detects  shocks. 
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I'or  ;i  inoi'c  driaiicl  i  iesrri  pt  k  m  ot  tlif  mult  i-dir<‘i’i  jonal  appioaiii  ilir  ii'adt'i'  ’i!a\ 
lie  ri'ten'ed  to  feiertMiee  in  these  proei'edin'i''^. 

Aftei'  a  deeade  of  miilti-iiireetioiial  methods,  wliat  iienelits  liav'e  heea  ileinoti- 
.•'tr.'ited?  Surely,  tliese  me'thods  yield  impressive  rt'siilts  u'Ik'u  ajydied  to  first-order 
seiiemes:  shock  and  shear  \va\'es  not  ali.u,iied  with  the  tirid  are  reprt'seiil ed  as  if 
coiitpiiti'd  with  a  hi!i'her-<vrd<u'  iriethod.  The  im[rrovei!ieiir  l(rom.iht  to  hiyiier-order 
schemes,  though,  is  a  lot  less  spectacular,  and  this  is  underst andahle.  On  tlu'  one 
hand,  there  is  not  much  room  ltd't  for  a  further  reduction  of  wa\e  spread  (mort>  for 
shtuir  waves  than  for  shock  waves  i;  on  the  other  hand,  loss  of  nionotonicity  may  occur, 
ayainst  which  there  are  no  etFecti\e  limiters,  ainl  convergenct'  to  a  steady  stale*  suffers 
umler  the  stiauig  nonlinearity  of  the  methods. 

In  m\’  opinion,  the  multi-diree't ional  approach  has  had  a  clear  itnjiae't  on  coniputa- 
rional  fluiif  dynamics.  .Although  complete  multi-directional  me'thods  will  sur\i\'e'  euily 
if  the  problem  of  ensuring  rejbustness  caii  be  solveel.  1  exjtee-t  that  ele'inents  of  siiedi 
me'thods  may  Knel  their  way  into  sfandarel.  elirection-split  codes,  to  help  re*solve  flow 
features  arising  in  specitie  iienv  prejl>le'ms. 


4.2  Minimum-strength  wave  models 

In  the'  work  of  Rumsr'y.  Reie  and  Van  Leer  [14]  anel  Parpia  and  .Michalek  [17].  the' 
orie'iitatiou  e)f  the  e-ell  interfae-e  is  de-emphasizeel.  The  spatial  discretization  is  no 
leuiger  regardeel  as  generating  a  dise-ontinuity  along  the  interfaces:  insteaei.  an  attempt 
is  made  to  find  out  what  waves  are  actually  propagating  near  the  interface.  This,  of 
e'ourse.  reejuire's  data  spanning  a  multi-dimensie)nal  ])art  of  space:  if  only  the  two 
states  I'i  and  I  r  are  to  be  used,  a  theoretieud  conjecture  must  make  up  for  the 
missing  information. 

In  the  basic  wave  model  of  Rumsey  et  al.  a  special  set  of  4  wavi’s  is  used  to 
match  the  state'  difference  I'r  —  I'l'.  for  uniqueness,  the  sum  of  the  wave  strengths 
is  miiiimiz('d.  Three  of  these  waves  follow  from  solving  a  one-dimensional  Riemann 
problem  in  the  direction  of  the  velocity  <lilference  Al  .  the  fourth  wave  is  a  shear 
wave  normal  to  the  other  three.  This  clioice  of  waves  makes  sense  from  a  kinematic 
[aunt  of  view,  as  illij.strat('d  by  Figure  10.  It  shows  that  a  velocity  difference  Al  can 
be  e.xplained  by  an  acoustic  wave  traveling  in  the  direction  of  A I  as  well  as  a  shear 
wa\e  traveling  in  the  normal  direction.  Which  exjtlanation  is  the  morr'  likeK'  one 
ma\'  be  determined  b\'  also  considering  the  |)ressure  difference  p/y  —  /qy  a  large  \alue 
favors  the  acoustic  explanation,  while  a  small  vahu'  favors  the  shear  ex|danation.  The 
minimization  procedure  takes  the  full  state  difference  (  h  —  I  l  and  comes  up  with  a 
[ilausible  ex])lanation  terms  of  all  four  waves.  The  method  of  Parpia  and  .Michalek 
differs  only  in  the  choice  of  the  functional  that  is  minimized.  Figure  II  shows  tin* 
configuration  of  the  plane  waves  crossing  the  interface.  In  practice  lurth  methods 
include  a  Hfth  wave,  a  weak  shear  wave,  whicli  corrects  for  the  difference  between  the 
true  direction  of  Al  and  the  direction  actually  used:  the  latter  may  have  been  held 
ov('r  from  a  previous  it  '»-ation  ("frozen"),  for  improvement  of  convergence. 

The  word  "plausible"  used  above  indicates  tiiat  the  minimization  procedure  only 
makes  an  educated  guess:  it  is  possible  to  com|)Ose  a  set  (jf  initial  values  that  is  totally 
misinterpreted.  ( 'onsiih'r.  for  insiam'e.  the  heail-on  collision  ol  two  gases  that  ha\(' 
e(]ual.  negligible  pressures.  In  rt-ality  two  strong  shocks  are  formed,  mo\  ing  into  t  lit' 


If) 


Figure  10:  Shock  or  shear? 

gases.  The  procedure  sees  as  input  a  velocity  difference  not  accompanied  by  a  pressure 
difference,  hence  calls  for  a  single  shear  wave,  as  if  the  gases  avoided  collision! 

The  flux  formula  based  on  the  above  wave  model  is  worth  some  discussion.  As¬ 
suming  the  system 

Ut  +  F{U),  +  G{iny  =  0,  (23) 

with  flux  .Jacobians  A{U)  and  represents  the  two-dimensional  Euler  equations, 

we  may  again  write  AT’  as  a  sum: 


SU  =  j2okRk.  (24) 

h=i 

The  vector  Rk  is  now  an  eigenvector  of  the  matrix 

A  cos  8k  4-  B  sin  Ok,  (25) 

where  9k  indicates  the  propagation  angle  of  the  A:-th  wave;  the  matrices  A  and  B  are 
standard  Roe-averages.  The  upwind-biased  interface  flux  is  defined  by 

1 

F{Ul,  Ur)  =  Fr)  -  ^  \\k  COii(0k  -  6^„om.al)|  f^kRk,  (26) 

^  t-1 

i.e.  still  by  formula  (21),  but  with  the  wave  speeds  A^,.  projected  onto  the  interface 
normal.  Although  this  formula  seems  trivial,  it  must  be  pointed  out  that  there  no 
longer  exists  a  relation  between  AF  and  AU  like  (6). 

In  numerical  practice  minimum-strength  wave  models  appear  to  bring  the  same 
benefits  and  problems  as  multi-directional  methods:  great  improvements  in  shock  and 
shear  resolution  for  first-order  methods,  much  smaller  improvements  for  second-order 
methods,  and  possible  loss  of  monotonicity  and  convergence. 

To  illustrate  the  performance  of  this  class  of  methods,  consider  Figures  Tia  and 
Fib.  Both  show  pressure  plots  for  steady  viscous  flow  over  ?  i\'.-\CA  0012  airfoil  at 
3°  angle  of  attack  and  Reynolds  number  5000,  computed  on  a  129  x  49  0-grid  by 
Rumsey  [12,  14].  Under  these  conditions  the  flow  separates  from  the  upper  surface. 
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Figure  11:  Plane  waves  crossing  a  cell  face  according  to  the  model  of  Rumsey  et  al. 

producing  a  detached  shear  layer  oblique  to  the  grid.  For  the  results  of  Figure  l‘2a  a 
second-order  MUSCL-type  scheme  [26.  44]  was  used,  with  Roe’s  [23]  standard  grid- 
aligned  Riemann  solver.  The  Riemann  solver  misinterprets  the  oblique  shear  as  an 
grid-aligned  shear  plus  an  acoustic  wave  (see  Figure  6);  the  latter  causes  a  pressure 
rise  or  drop  at  tlie  interface.  Correspondingly,  the  steady  solution  shows  pressure 
fluctuations  across  the  shear  layer,  so  that  its  presence  can  actually  be  detected  in 
pressure  plots.  .A  grid-refinement  study  shows  that  the  disturbances  scale  with  the 
mesh  size.  This  phenomenon  was  first  observed  by  Venkatakrishnan  [41]  and  correctly 
explained  by  Roe:  in  fact,  it  motivated  the  work  of  Rumsey,  Van  Leer  and  Roe.  As 
seen  from  Figure  12b.  the  minimum-strength  wave  model  properly  recognizes  the 
oblique  shear  layer  and  generates  clean  pressure  contours. 

The  same  method  gives  an  unexpected  improvement  in  the  representation  of  in- 
viscid  stagnating  flow.  The  explanation  is  found  in  Figure  13,  showing  the  turning 
of  the  flow  near  a  stagnation  point  S  cis  represented  by  the  discrete  velocities  in  the 
three  cells  marked  1,  2  and  3..  A  grid-aligned  Riemann  solver  interprets  the  velocity 
ditfevence  between  vertical  neighbors  1  and  2  as  a  compression  (V^i  >  Vy2),  and  the 
velocity  difference  between  horizontal  neighbors  2  and  3  as  an  expansion  (142  <  ^43); 
this  leads  to  pressure  variations  of  the  order  of  AV.  The  wave  model  detects  only 
very  small  pressure  changes  (Ap  ~  pA(V'^))  and  therefore  explains  both  velocity  dif¬ 
ferences  by  shear  waves.  Although  this  still  is  not  the  right  explanation,  the  result  is 
a  decrease  in  numerical  entropy  production.  The  effect  is  rather  large  for  first-order 
methods,  as  can  be  judged  from  Figure  14  showing  entropy  contours  for  inviscid  flow 
over  a  NACA  0012  airfoil  at  M  =  0.3,  a  =  1°,  on  a  sequence  of  0-grids.  The  reduced 
entropy  levels  lead  directly  to  reduced  numerical  drag  levels,  eis  Figure  15a  shows.  For 
second-order  schemes  the  effect,  cis  usual,  is  less  dramatic;  the  drag  values  are  given 
in  Figure  15b. 
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Figure  12:  Viscous  separating  flow  over  a  NACA  0012  airfoil  at  M  =  0.5,  a  =  3°  and 
Re  =  5000.  Pressure  contours  on  a  129  x  49  C-grid,  obtained  with  a  second-order 
upwind  scheme  incorporating  (a)  Roe’s  grid-aligned  Riemann  solver;  (b)  the  five-wave 
model  of  Rumsey  et  al. 


5  Multi-dimensional  fluctuation  approach 


The  fluctuation  approach  to  upwind  differencing  lends  itself  better  to  e.xtension  into 
higher  dimensions  than  the  finite-volume  approach.  Recall  that  a  fluctuation  is  a 
local  flux  imbalance  causing  a  non-zero  time  derivative  of  the  local  solution.  For  the 
one-dimensional  Euler  equations  (4)  the  quantity  —AF  equals  the  residual  evaluated 
on  a  one-dimensional  mesh: 


I  Utdx  =  -  I  F^dx  =  -AF. 

J  mesh  J  mesh 


(27) 


This  suggests  extension  of  the  fluctuation  approach  beyond  one  dimension  by  regard¬ 
ing  each  multi-dimensional  mesh  residual  as  the  sum  of  a  finite  number  of  vv^aves  (say, 
ni),  moving  in  all  possible  directions.  Thus  we  discretize  the  two-dimensional  Euler 
equations  a.s 


/  /  Utdxdxj  =  -  fiFdij  -  Gdx)  -  A^aA,./?*,., 

J  Viuesh  J  k=\ 


(28) 
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Figure  13:  Turning  of  the  flow  in  three  cells  near  a  stagnation  point  5  at  a  wall. 


(a)  (b) 


Figure  14:  Entropy  contours  for  inviscid  flow  over  a  NACA  0012  airfoil  at  M  =  0.3, 
a  =  1°.  generated  on  a  sequence  of  0-grids  with  a  first-order  scheme  incorporating 
(a)  Roe's  grid-aligned  Riemann  solver;  (b)  the  five-wave  model  of  Rumsey  et  al. 

where  the  matrices  A  and  B  are  multi-dimensional  averages  that  remain  to  be  de¬ 
fined.  .Since  the  fluctuation  approach  is  a  nodal-point  approach,  and  we  wish  to 
develop  only  schemes  of  maximum  compactness,  we  shall  use  a  grid  of  triangular 
meshes,  with  data  given  in  the  nodal  points.  For  the  computation  of  the  residual 
on  such  meshes  it  suffices  to  apply  the  trapezoidal  integration  rule  on  each  side  of 
the  triangle.  The  fluctuations  resulting  from  residual  decomposition  must  be  sent  to 
the  triangle's  vertices  according  to  some  distribution  scheme  that  approximates  the 
convection  equation. 

It  follows  that,  for  the  construction  of  a  genuinely  multi-dimensional  upwind- 
differencing  scheme,  three  components  are  needed: 

1.  A  reliable  multi-dimensional  wave  model  for  representing  the  residual; 

2.  .A  way  to  ensure  conservation,  i.e.  a  multi-dimensional  extension  of  Roe’s  matrix 
average: 


1/(n:-nj)’^  '/(ni*nj) 


Figure  15:  Grid-convergence  study  of  the  drag  coefficient  based  on  (a)  the  first-order 
solutions  of  Figure  14;  (b)  the  corresponding  second-order  solutions. 

3.  .A.  multi-dimensional  convection  scheme  for  advancing  the  waves. 

Each  of  these  will  be  discussed  in  a  separate  subsection. 

5.1  Multi-dimensional  wave  models 

The  modeling  of  a  local  Euler  residual  by  a  finite  number  of  waves  was  launched  as 
a  research  subject  by  Roe  [2];  his  first  paper,  however,  gave  no  specific  instructions 
as  to  how  the  model  would  be  used  in  a  numerical  integration  of  the  Euler  equations. 
This  is  not  surprising,  given  that  the  other  problems  -  multi-dimensional  conservation 
and  advection  -  had  not  yet  been  addressed. 

The  latest  version  of  Roe's  wave  model  calls  for  four  acoustic  waves,  running  along 
the  principal  strain  axes  of  the  local  fluid  element,  a  shear  wave  making  a  45®  angle 
with  the  acoustic  waves,  and  an  entropy  wave  running  in  the  direction  of  the  entropy 
gradient;  see  Figure  16.  Thus,  m  =  6  in  Eq.  (28).  These  six  waves  are  defined 
by  two  independent  angles  and  six  strengths:  therefore,  eight  independent  pieces  of 
information  need  to  be  supplied  per  triangular  mesh.  This  information  is  available 
in  the  form  of  the  gradient  of  the  state  vector;  its  mesh  value  is  computed  with  the 
trapezoidal  rule  from  the  following  boundary  integrals: 

f  X  =  -T^ —  f  f  b\dxdy  =  -j^  /  Udy.  (29) 

AlJ'CCI  J  ^niesli  t(l  ^mesh 

Uy  =  —  /  /  Uydxdy  =  — I  Udx.  (30) 

AvCd  J  Vmesli  ArCd  2niesh 

A  detailed  discussion  of  this  wave  model,  including  the  three-dimensional  case,  can 
be  found  in  Roe’s  contribucion  to  the  present  volume  [45];  numerical  results  obtained 
with  this  model  are  presented  in  the  contribution  by  Catalano  et  al.  [46]. 

This  section  would  not  be  complete  without  a  discussion  of  the  work  of  Hirsch  and 
collaborators  [47,  48,  49].  Their  multi-dimensional  approach  is  based  on  diagonalizing 
the  Euler  equations,  i.e.  changing  these  into  a  system  of  convection  equations,  by  a 
transformation  of  state  variables.  The  transformation  itself  depends  on  the  local 
gradient  of  the  solution,  making  the  diagonalization  essentially  nonlinear.  For  certain 
data  the  transformation  does  not  exist,  in  which  case  it  is  chosen  so  as  to  minimize 
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Figure  IG:  Roe  s  two-dimensional  six-wave  model.  The  acoustic  waves  run  parallel 
to  the  principal  strain  axes  (dashed);  the  strain  ellips  (dotted)  shows  the  kinematic 
deformation  of  a  circular  fluid  element. 

the  otf-diagonal  terms.  The  update  scheme,  though,  can  be  made  identical  to  a 
fluctuation-based  scheme:  decomposition  of  the  residual  along  certain  eigenvectors, 
followed  by  convection  of  the  components  [oO].  In  two  dimensions  the  diagonalization 
is  equivalent  to  using  one  particular  four-wave  model:  clearly,  the  fluctuation  approach 
offers  much  more  flexibility. 


5,2  Multi-dimensional  conservation 

The  multi-dimensional  extension  of  Roe's  averaging  of  the  flux  .Jacobian  was  indepen¬ 
dently  discovered  by  Roe  and  Struijs.  and  is  presented  in  a  joint  paper  [51].  This  very 
recent  (1991)  addition  to  the  multi-dimensional  toolbox  applies  exclusively  to  trian¬ 
gular  meshes  in  two  dimensions  and  tetrahedral  meshes  in  three  dimensions.  The 
following  (.lescription  and  explanation  of  the  two-dimensional  averaging  apply  to  the 
special  case  of  a  calorically  perfect  gas. 

To  begin  with,  assume  that  the  parameter  vector  w  =  v.  H]  is  distibuted 

linearly  over  a  mesh  triangle  with  vertices  labeled  1.  2  and  3.  Denote  the  average  of 
u'  over  the  triangle  by  il".  we  then  have 

W  —  -(U’l  -h  ti'2  -t-  U’3).  (31) 

.As  before,  l'(w)  and  F{nj).  and  also  are  quadratic  in  the  components  of  il\  so 

that  the  .Jacobian  matrices  (  „,.  and  Gu,  are  linear  in  vo.  and  therefore  also  in  x 

and  y.  Considering  that  I'j.  =  I'y  =  etc.,  where  u;,.  and  Wy  are  constant 

(jver  the  entire  triangle,  we  conclude  that  VC,  VC  and  VG'  also  vary  linearly  over  the 
triangle,  ('sing  the  definition  of  the  mesh-averaged  gradient  Vf'  given  in  Eqs.  (29), 
(30),  and  similar  definitions  of  VF  and  VG',  we  easily  derive  the  relations. 


Figure  17;  Stencils  of  two-dimensional  upwind  convection  schemes;  case  n  >  6  >  0. 

(a)  Sidilkover's  second-order  scheme.  The  fluxes  for  cell  1  nominally  are  compurted 
by  linear  interpolation  between  upstream  pairs  of  data,  but  the  fluxes  at  the  North 
and  South  faces  must  be  limited  to  prevent  numerical  oscillations.  The  limiters  are 
based  on  the  ratios  a(ui  —  U2)/[b{us  —  U2)]  and  a{u^  —  ■U4)/[6(u2  —  U4)],  respectively. 

(b)  Standard  second-  or  third-order  grid-aligned  scheme. 


which  are  direct  extensions  of  the  one-dimensional  relation  (6).  The  extension  to 
three-dimensional  averaging  is  self-evident. 


5.3  Multi-dimensional  convection 

The  pursuit  of  multi-dimensional  convection  schemes  has  kept  a  number  of  authors 
busy  over  the  past  three  years.  In  two  dimensions  the  basic  equation  to  be  solved  is 

Ui  -I-  aui  +  buy  =  0.  (34) 

where  a  and  b  are  constant  velocity  components,  or,  in  vector  notation. 

u;  4-  a  •  Vu  =  0  (35) 

The  first  significant  work  was  that  of  Sidilkover  [52],  who,  among  other  things, 
showed  how  a  second-order  upwind  scheme,  with  residual  computed  on  a  square  mesh, 
can  be  made  non-oscillatory  by  standard  limiters  without  undue  spreading  of  the 
stencil.  The  domain  of  dependence  for  this  algorithm  is  shown  in  Figure  17a.  for  the 
case  a  >  b  >  0\  note  how  compact  this  is  in  comparison  to  the  stencil  of  a  standard 
second-order  upwind  scheme,  shown  in  Figure  17b  [27].  He  also  coined  the  name 
”N-scheme”  for  the  first-order  scheme  that,  on  a  cartesian  grid,  takes  its  data  from 
the  upwind  triangle  fitting  the  convection  path  most  tightly  (N  stands  for  narrow). 
For  example,  for  point  1  in  Figure  17a  it  would  be  triangle  (124).  This  scheme,  as 
shown  in  [53],  is  optimal  in  the  sense  that,  among  all  schemes  with  upwind  triangular 
domain  of  dependence,  it  combines  the  smallest  truncation  error  with  the  largest 
stable  time-step.  The  three-dimensional  extension  is  also  described  in  [53]. 

While  the  triangles  in  Sidilkover’s  work  were  still  considered  subdivisions  of  squares, 
they  become  autonomous  in  later  work  by  other  authors.  A  major  step  in  the  devel¬ 
opment  of  two-dimensional  convection  schemes  was  the  realization  that  there  are  two 
types  of  triangles  [54]:  those  with  one  inflow  side  and  those  with  two  inflow  sides. 
This  is  illustrated  in  Figure  18.  If  there  is  only  one  inflow  side,  the  fluctuation  aj)- 
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i'iiiurt'  IS:  Two  kinds  ot  trian“;les:  (a)  with  one  intiow  si(ie:  (b)  with  two  iiiHow  sides. 

proai'h  dictates  that  the  entire  residual  be  used  to  update  the  opposite  node.  This  is 
the  unique  ".single-tar2;et"  t'orni  of  the  scheme,  similar  to  the  oue-dimensiouai  upwind 
s<  heme  2.  If.  however,  tliere  are  two  iuHow  sides,  it  may  be  argued  that  the  residvial 
be  distril)utt'd  over  the  two  nodal  points  defining  the  third  side.  This  is  the  "dual¬ 
target"  form  of  the  particular  scheme;  each  choice  of  distribution  weights  defines  a 
new  si'heme.  The  s})reading  of  the  residual  information  over  two  points  implies  a 
potential  loss  of  r*  solution,  inlierent  to  multi-dimensional  numerical  convection:  there 
is  no  one-dimensional  analogU(»  of  this  elfect. 

In  the  develoijment  of  nndti-dimensional  convection  schemes,  three  design  criteria 
play  a  decisive  role,  .-\ccording  to  these,  it  is  desirable  for  a  scheme  to  be 

I .  linear:  for  a  given  grid  geometry  ainl  flow  angle  the  solution  depends  linearly  on 
the  flat  a.  This  prontotes  <'onvt>rg<'nce  to  a  steady  numerical  solution.  It  is  well 
km)wn  that  tht'  [)resence  of  nonlinear  devices  in  the  scheme,  such  as  limiters  [4-  j 
and  frame  rotation  (see  .Section  4.2)  can  slow  down  or  even  halt  the  convergence 
process: 

J.  linearity  preserving  (LP):  data  of  the  form 

u(.r,  //)  =  hx  —  (ly.  (db) 

which  is  a  steady  solution  of  Eq.  44  are  not  changed  by  the  scheme.  This 
promotes  the  accuracy  of  tin’  s('heme.  It  can  be  shown  [54]  that  LP  schemes 
yi '  Id  second-order-accurate  steady  s<.»lutions  ol  E<|.  -M: 

.5.  positive:  the  scheme  has  positive  coefficients.  This  is  sufficient  for  preventing 
numerical  (jsci  Hat  ions. 

From  one-dimensional  finite-difference  theory  we  know  -  and  have  known  so  for 
a  long  time  -  that  the  abme  conditions  are  mutually  exclusive.  There  is  a  famous 
theoieni  by  (lodunov  [24]  which  says  that  no  litiear  convection-diffusion  scheme  with 
positi\e  coefficients  can  be  more  than  first-order  accurate.  With  reference  to  our  de¬ 
sign  criteria  for  multi-dimensional  cotivecfion  schemes  this  theorem  reads: 

Thrvf  an  an  liat  arpnsitirr  L  f*  srhi  ait  s. 

.Again,  noidinearity  is  essential  for  the  design  of  accurate,  non-oscillatory  schemes. 


IS 
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(a)  (b) 


Figure  U);  Dual-target  form  of  convection  schemes;  (a)  :V-scheme;  (b)  LDA-scheme, 

Among  the  various  upwind  convection  schemes  proposed  in  recent  years,  three 
s('hemes  stand  out:  these  are  discussed  l)elow.  They  all  are  as  compact  as  can  be. 
requiring  data  on  only  one  triangle  for  the  approximation  of  the  convection  equation. 
A  small  miracle  is  that  even  positivity  can  be  a<'hieveil  without  having  tht  trianglt. 
Of  course,  each  norlal  point  is  a  vertex  of  a  number  of  triangles  and  may  receive 
fluctuations  from  several  ot  these;  programming  therefore  must  be  triangle-based. 
Some  results  of  numerical  experiments  are  presented  in  Section  6. 

The  N-scheme:  the  optimal  linear  positive  scheme 

The  name  of  this  scheme  suggests  equivalence  to  Sidilkover's  .\'-scheme,  but  it  actu¬ 
ally  is  more  general.  .Sidilkoverbs  scheme  is  just  the  single-target  form,  common  to 
all  compact  schemes:  Huctuations  from  triangles  requiring  a  dual-target  scheme  are 
ignored  in  the  up<lale.  The  dual-target  form  of  the  current  .\-scheme  uses  distribution 
weights  proimrtional  to  the  componetits  of  the  convection  speed  along  the  two  inflow 
sides,  as  depicted  in  Figure  19a.  This  makes  the  scheme  optimal  in  the  sense  of  having 
the  largest  stability  range  for  the  time-step  [-54].  It  is  also  linear  and  positive,  and 
therefore  can  be  no  more  than  first-order  accurate. 


The  NN-scheme:  the  optimal  nonlinear  positive  LP  scheme 

This  scheme  is  a  nonlinear  variant  of  the  .\-scheme.  hence  the  second  .N' .  The  nonlinear 
procedure  included  in  this  scheme  has  alxsointely  nothing  in  common  with  the  TV'D- 
enforcing  limiters  included  in  one-dimensional  convection  schemes.  It  is  based  on 
the  observation  that  in  the  convection  ecpiation  (35)  the  component  of  the  convection 
velocity  a  perpendicular  to  the  solution  gradient  V«.  has  no  effect  on  U(.  We  therefore 
are  allowed  to  replace  a  by  any  velocity  that  has  the  same  component  parallel  to  Vu. 
as  shown  in  Figure  20.  This  compotient.  indicated  by  n,^..  is  the  velocity  at  which 
the  level  lines  of  u  normal  normal  to  themselves,  i.e.  the  wave  speed  of  the  local 
distribution  of  u.  This  wave  speed  is  the  smallest  of  all  admissible  convection  speeds: 
it  actually  vani‘=hcs  with  the  residual.  We  may  now  adopt  the  following  strategy;  if 
both  a  and  a„.  call  for  a  (hial-target  sclumie.  we  reiilace  u  by  (7„,  in  the  N-scheme;  in 
all  other  cases  the  schtmie  becomes  or  remains  a  siiigh'-target  scheme.  In  the  case  t)f 


If) 


Figure  20;  The  .\’.\-.srhenie:  nonlinear  single-target  form.  The  convection  speed  a 
calls  for  a  dual-target  scheme  but  is  replaced  by  aj,  which  calls  for  a  single-target 
scheme.  The  wave  speed  d,u  is  the  component  of  a  and  d]  parallel  to  Vu. 

Figure  20.  a  is  replaced  by  Si.  the  nearest  admissible  speed  yielding  a  single-target 
scheme.  The  resulting  scheme  does  not  change  any  nodal  value  if  the  residual  vanishes, 
hence  is  LP,  and  maxindzes  the  allowable  time  step. 

.N’umerical  results  indicate  that  the  accuracy  of  the  NN-scheme  lies  between  first- 
aiid  second-order:  see  further  .Section  6. 


The  LDA  scheme:  a  non-positive  linear  LP  scheme 

This  scheme  is  one  of  several  low-diffusion  schemes,  designed  for  a  low  truncation 
error.  In  tlie  dual-target  form  of  the  scheme  the  distribution  weights  are  inversely 
proportional  to  the  areas  of  tlie  triangles  cut  from  the  mesh  triangle  by  a  streamline 
through  the  inflow  vertex:  see  Figure  19b.  This  scheme  is  not  positive,  but  very  accu¬ 
rate:  on  a  uniform  gritl  it  achieves  third  order  accuracy,  as  demonstrated  in  Section 
(). 


The  above  schemes  have  served  as  the  basis  for  convection-diffusion  schemes  in  a 
study  l)y  Tomaich  and  Roe  [o')].  .Since  the  diffusion  operator  can  not  be  approximated 
on  a  single  triangle,  their  schemes  are  formulated  with  reference  to  a  central  nodal 
point.  Numerical  solutions  of  the  Smith-Hutton  [-oG]  test  problem  demonstrate  that 
these  schemes  rival  the  best  exsisting  convection-diffusion  schemes  in  accuracy.  In 
a<ldition.  their  way  of  discretizing  the  Laplacean  is  directly  applicable  to  any  of  the 
disipative  terms  included  in  the  Navier-.Stokes  equations;  thus,  the  basis  for  genuinely 
midti-dimensional  N'avier-Stokes  codes  has  been  laid. 


6  Numerical  results 

To  support  some  of  the  statements  made  about  the  new,  compact  convection  schemes 
I  first  show  how  these  schemes  fared  in  a  comparative  grid-refinement  study  by  .Jens 
.\Iiiller  [.')7|.  The  problem  is  that  of  convection  of  a  Claussian  distribution  over  a  semi¬ 
circle;  Inflow  is  at  //  =  0.  .r  <  0.  outfluw  at  (/  =  0.  x  >  0.  Four  kinds  of  grids  were 
used,  of  which  three  examples  are  displayetl  in  Figure  21.  Grids  a  and  3  derive  from 
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ri‘iUM'  l!  1 :  1  liirt*  lirius  list’d  iii  the  rireular-cenvi'i't ion  ex])cnnieiits. 


v-3 - --^3  N  oil  <1 

— c  N  oil  J 

3  •  -  -  c  N  oil  7 
13-  •  —  23  N  on  <5 

<5 - «  NN  Of)  a 

NN  on  J 
®  -  -  -  •  o  N  N  oil  7 
o - o  NN  on  6 

6 - A  LHA  on  o 

ft - ft  LDA  on  J 

ft - ft  LDA  on  7 

ft—  •  —  ft  L  D  A  on  6 
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I'  inure  JL’;  ( iri<i-<'(j[i\  <’r2;(’nr(’  of  Z,j-err<>r  tor  e(;n\'e('i ion  of  a  (  jau-s.^ian  over  a  .senii('ir('le 
liy  \'arious  scluMiies  on  \arious 

a  unitorni  cartesian  nrid  h\’  aildinir  diaiionals.  in  'zrid  f>  those  dianonals  are  ehosiui 
that  are  h'ast  aliniied  with  tht'  euiiveetion  direetmn.  in  J  those  most  alietied.  (iiid 
'  is  a  irregular  perturhatiuti  to  d.  wliile  d  (not  shown)  is  a  minor  jun't  m  hat  ion  to 
■.  Fisure  dll  shows  the  eonv(U'2;en<<’  ol  the  /ij-error  produced  hy  t  lie  X-.  XX-  and 
LDA-sehemes  on  the  rhtFerent  grids.  From  tlie  slope  of  the  graphs  of  |og(errorl  versus 
lon(  mesli- width )  the  following  roindiisions  can  lie  drawn; 

I.  Fhe  ,X-seheme  is  som<nvhat  h’ss  than  first-order  accurate; 

d.  Fhe  X.X-scheme  is  cdoser  to  hi’ing  secoml-order  ai'ciirate  than  tirst-(,)rder  accu¬ 
rate; 

d.  Fhe  Ll)A-s(  heme  is  third-order  accurate  on  a  renuiar  yrid.  second-order  or  les;' 
on  a  perturlied  grid; 


'ft 


Fissure  'J-F  Contours  (left)  and  cut  along  the  a‘-axis  (right)  of  the  solution  obtained 
with  the  LDA-srheine  on  a  20  x  20  .Cgrid. 


Figure  24:  Iriviscid  How  at  Mach  1.4  through  an  inlet,  computed  by  a  genuinely 
multi-(limensional  upwind  Euler  code.  .Shown  are  Mach-nuinber  contours  and  the 
unstructured  grid  used. 

4.  .All  schemes  decrease  their  error  when  diagonals  are  aligned  with  the  How. 

.Most  suri)rising  is  the  achievement  of  tliir<l-order  accuracy  on  regular  grids,  consider¬ 
ing  the  limited  amount  of  information  going  into  these  compact  scliemes.  Figure  23 
gives  an  idea  of  this  high  accuracy  by  showing  solution  contours  and  a  cut  at  ;/  =  0 
obtained  with  the  LD.A-scheme  on  the  very  coarse  /i-grid  of  Figure  21  (Gaussian  rep¬ 
resented  on  10  meshes).  .Similar  results  obtained  with  three-dimensional  extensions 
of  the  schemes  ran  be  found  in  Deconinck's  comprehensive  review  paper  [4]. 

While  the  search  for  compact  convection  schemes  continues,  several  authors  are 
trying  to  put  together  the  ingredients  listetl  in  .Section  5,  producing  a  genuinely  multi¬ 
dimensional  upwind  Euler  code.  .Advanced  numerical  results  can  be  found  in  the 
present  volume  in  the  contribution  by  Catalano  et  al.  .An  earlier  successful  calcvda- 
ti(ui  (d  su[)ersonic  iidet  How  Iry  Struijs  et  al.  [3]  produced  the  Marli-contours  shown 
in  Figure  24:  superimposed  is  the  moderately  irregular  triangidation.  The  results 


Figure  25;  The  finite-volume  scheme  of  Barth.  Powell  and  Parpia.  based  on  a  trian¬ 
gular  grid  and  its  median  dual.  The  flux  F]  across  one  median  element  located  in 
triangle  T  is  computed  using  wave  data  from  triangle  T  and  fluxes  from  the  vertices 
Li  and  R.  For  the  flux  F^  across  the  other  element  in  the  same  triangle  the  same 
wave  data  are  used,  but  the  fluxes  are  taken  from  L2  and  R. 

demonstrate  the  excellent  shock-capturing  ability  of  the  NN-scheme. 


7  Finite- volume  schemes  revisited 

From  Section  5  the  reader  might  get  the  impression  that  genuinely  multi-dimensional 
schemes  can  only  be  formulated  on  triangular  and  tetrahedral  meshes,  and  that  they 
are  incompatible  with  the  finite-volume  formulation.  If  this  were  true,  it  would  mean 
a  serious  restriction  on  the  use  of  such  schemes  within  the  CFD  community,  for  it 
is  not  at  all  clear  that  unstructured  triangular  or  tetrahedral  grids  are  the  way  of 
the  future.  An  alternative,  for  instance,  is  offered  by  adaptive  cartesian  grids  [58]. 
Tlie  emphasis  on  triangles  in  Section  5  arises  from  the  experience  that  the  numerical 
building  blocks,  e.g.  Roe's  matrix  averaging,  take  their  simplest  form  on  such  meshes. 
Therefore,  in  developing  a  multi-dimensional  scheme  of  a  different  format  it  would  be 
good  practice  to  start  with  the  wave  decomposition  of  residuals  on  an  underlying 
triangular  grid. 

As  an  example  of  such  practice,  consider  the  two-dimensional  finite- volume  Euler 
scheme  of  Powell.  Barth  and  Parpia  [59],  illustrated  in  Figure  25.  The  wave  model 
indeed  is  applied  to  data  on  triangular  meshes;  for  the  update,  though,  a  finite- volume 
scheme  is  chosen.  Cell  faces  are  formed  from  the  medians  of  the  triangles,  yielding 
the  so-called  median-dual  grid.  Across  each  median  element  of  the  cell  contour  a  flux 
is  computed  using  an  equation  of  the  form  (26),  where  L  and  R  denote  the  vertices 
of  the  triangle  side  bisected  by  the  median,  and  the  sum  includes  all  waves  identified 
in  the  triangle.  Note  the  difference  with  the  scheme  of  Rumsey  et  ah,  where  the  wave 
model  would  be  ba.sed  solely  on  Ur  —  iU-  The  resulting  nonlinear  scheme,  applied 
to  a  scalar  convection  equation,  is  LP  and  positive  and  appears  to  be  more  accurate 
than  the  NN-scheme. 
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8  Conclusions 


I  !ic  siatf  ol  tlu-  art  in  ,y;enuinely  nmlti-rliinonsional  upvvin<1  differrurinir  has  made 
di';imatic'  aiKatua's  over  the  past  throt'  years,  owin*;;  tu  a  shift  from  the  finit(’-\’oliim(‘ 
appro.u'ii  to  the  lict nation  approach.  The  basic  ingredients  for  inulti-dirntnisiunal 
I'.nlcr  codes,  i.e.  wave  model,  conservation  principle  and  convi'ction  scheme,  are  read\- 
lor  intt-gration.  and  f.he  first  numerical  results  look  good.  The  coming  years  will  viidd 
maiiv  more  Eulei'  applicatiiuis  in  two  and  three  dimensions,  further  improvements  in 
wa\e  motlels  and  compact  convection  schemes,  and  extension  of  the  approach  to  tlu' 
iiujdeling  of  the  .Xavier-Stokes  equations. 
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